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Abstract 

In this work we show how the inclusion of dark matter (DM) direct detection upper 
bounds in a theoretically consistent manner can affect the allowed parameter space of a 
DM model. Traditionally, the limits from DM direct detection experiments on the elastic 
scattering cross section of DM particles as a function of their mass are extracted under 
simplifying assumptions. Relaxing the assumptions related to the DM particle nature, 
such as the neutron to proton ratio of the interactions, or the possibility of having similar 
contributions from the spin independent (SI) and spin dependent (SD) interactions can 
vary significantly the upper limits. Furthermore, it is known that astrophysical and nu¬ 
clear uncertainties can also affect the upper bounds. To exemplify the impact of properly 
including all these factors, we have analysed two well motivated and popular DM sce¬ 
narios: neutralinos in the NMSSM and a Z' portal with Dirac DM. We have found that 
the allowed parameter space of these models is subject to important variations when one 
includes both the SI and SD interactions at the same time, realistic neutron to proton ra¬ 
tios, as well as using different self-consistent speed distributions corresponding to popular 
DM halo density profiles, and distinct SD structure functions. Finally, we provide all the 
necessary information to include the upper bounds of SuperCDMS and LUX taking into 
account all these subtleties in the investigation of any particle physics model. The data 
for each experiment and example codes are available at this site http: //goo. gl/lCDFYi, 
and their use is detailed in the appendices of this work. 


1 Introduction 


Nowadays, the dark matter (DM) nature is one of the greatest mysteries of modern physics. 
The possibility of non gravitational interactions of DM with ordinary matter is one of the key 
pieces to understand the DM problem. In this sense, the experimental community is putting 
a great deal of effort to provide such a valuable information. 

Direct searches for DM aim at detecting the scattering of DM particles off nuclei inside 
a low-background target detector. Once a DM particle hits a nucleus, the recoiling nucleus 
will release a certain amount of energy that, above a given threshold, can be measured using 
different and sophisticated techniques. This field is currently undergoing an exciting period, 
with several experimental collaborations, using different targets and techniques, reporting po¬ 
tential signals of DM that are, nevertheless, in conflict with null results in other experiments. 
Most remarkably, the results obtained by the SuperCDMS [1] and LUX [2] Collaborations 
have placed the most stringent upper bounds on the elastic scattering spin independent (SI) 
cross section up to date, for DM masses above 4 GeV. Particle physics models for DM are 
then tested using these results, which in general, reduce considerably the allowed region of 
the parameter space of a given model. 

We are living in the precision data era. There is a huge amount of experimental data 
from colliders, direct detection and indirect detection experiments, that must be taken into 
account to determine the viability of a model of physics beyond the Standard Model (BSM). 
In this regard, theoreticians are making a special effort to reduce the uncertainties from the 
theoretical side, for instance, in the calculation of the Higgs mass, which improves notably 
the interpretation of the experimental data in the context of different BSM scenarios. 

Direct detection experimental collaborations adopt a set of assumptions about the DM 
nature and the DM halo which are useful in order to compare different results within a 
unified framework. More specifically, their results are traditionally presented considering 
only SI interactions^ and the same interaction strength of DM with protons and neutrons, 
i.e. with a neutron to proton ratio fn/fp = 1- The speed distribution of DM particles is 
taken to follow a Maxwell-Boltzmann distribution, a consequence of assuming the halo to 
be an isothermal and isotropic sphere, the Standard Halo Model (SHM). Then, these results 
are used to constrain the parameter space of models, like supersymmetry (SUSY), in which 
theoretical calculations are performed with a high degree of accuracy. Furthermore, in many 
models the assumptions that experimental collaborations consider do not apply neither for 
the DM particle considered nor for the DM halo profile^. This is a consequence of the 

^Sometimes, they also present bounds considering only spin dependent (SD) interactions under certain 
assumptions about the neutron to proton ratio an jap. 

^For DM candidates with a mass heavier than approximately 100 GeV the impact of speed distribution 
uncertainties on the upper limits of the elastic scattering cross section can be generally neglected. 
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aforementioned over-simplifying assumptions, as realistic particle physics models of DM (e.g. 
SUSY) allow for different couplings to neutrons and protons and similar contributions from SI 
and SD interactions. Then, in principle, it is not fully consistent to employ the upper limits 
quoted by experimental collaborations to constrain general particle physics models of DM. 
However, to take into account all these subtleties is not an easy task, since it would require 
the simulation of a direct detection experiment, which makes very difficult their inclusion 
when scanning the parameter space of a model. 

All in all, it is usual to compare the predictions of a given model with direct detection 
upper limits at face value, this is, considering the SI and SD interactions separately^, an equal 
contribution from protons and neutrons, i.e. fn/fp = 1, and assuming the SHM. Nonetheless, 
it has been pointed out that either the neutron to proton ratios [3-8] or the SD structure 
functions [9] and speed distributions [10-16] chosen can vary the interpretation of direct 
detection data substantially. In fact, to account for the latter effect, several halo independent 
methods have been proposed recently [17-28]. In this article, we present a set of tabulated 
data that contain the necessary information to calculate the expected number of DM induced 
events in the SuperCDMS and LUX experiments taking into account different possibilities 
of the mentioned ingredients. These data allow a fast check of the validity of a given model 
solution considering both, the SI and SD interactions, for arbitrary values of fn/fp and an/ap, 
using realistic DM speed distributions and including the uncertainties related to the SD form 
factors. We provide the data and example codes at this site http://goo.gl/lCDFYi. 

To this end, we have carried out a decomposition of the SuperCDMS and LUX detector 
signals in intensity (the part related to the cross sections) and shape (the part dependent on 
the energy) in such a way that the total number of predicted events in these experiments 
can be easily calculated even if the particle physics and/or DM profile assumptions made 
by collaborations do not allow a consistent and serious test. The tabulated data we have 
derived are very similar to the scaling functions introduced in Ref. [29]^. Here, we extend 
this information to include more realistic descriptions of the DM halo profile, and different SD 
structure functions. These data can be used to evaluate at which confidence level a certain 
point of the parameter space of a model is excluded or allowed, without any computationally 
intensive calculation. Such an observable can then be translated into upper limits on the 
properties of the DM candidate. The aim behind this work is to provide the particle physics 
community with enough information to include direct detection bounds in terms of the DM 
expected number of events, in a complete and consistent manner. 

To exemplify the impact that the combined limits (using the SI and SD interactions at 
the same time), together with distinct assumptions about the speed distribution and the SD 

^Sometimes, not even considering the SD interactions for a spin 1/2 DM particle. 

■^Notice that the definition of these functions in Ref. [29] varies slightly respect to the case presented here. 
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structure functions, have on the allowed solutions in a model parameter space, we have chosen 
as case studies two well motivated and popular models: the Next-to-Minimal Supersymmetric 
Standard Model (NMSSM) and a Z' portal with Dirac DM. We show explicitly step by step 
how the previous considerations can affect the phenomenological viability of the solutions 
found. Moreover, the data provided in this work for different DM halo profiles, can be used 
to perform a consistent analysis of a model in light of DM direct and indirect detection 
experiments using the same DM halo density prohle assumptions. 

This paper is organised as follows. In Section 2, we present the basic expressions that 
allow to compute a DM signal in the LUX and SuperCDMS detectors. Then, in Section 3 we 
derive the upper limits from these experiments under different assumptions in the velocity 
distribution and SD form factors, emphasizing their impact on the predicted number of events. 
In Section 4 and 5, we show how the inclusion of the upper bounds in terms of the number 
of events, considering both the SI and SD interactions, realistic neutron to proton ratios and 
different speed distributions, has a notable impact on the allowed parameter space of the 
NMSSM and Z' portal models, respectively. Then, we repeat this exercise using different SD 
structure functions, in such a way that we can quantify the effect on the parameter space of 
these scenarios. Finally, the conclusions are left for Section 6. In appendix A, we give detailed 
information about the extraction of the upper bounds and how to apply them to a generic 
model, and in appendix B we briefly describe the data files and example codes attached to 
this work. 


2 Direct detection of dark matter 


Let us start by summarising the basic expressions that describe the differential rate in a DM 
direct detection experiment [30] (for a recent review, see Ref. [31]). The differential event 
rate for the elastic scattering of a WIMP with mass off a nucleus with mass mjsi reads® 

do 


dR 


Po 


dEji rriN rriy. 


vf{v) 


dEji 


[v, Er) dv , 


( 2 . 1 ) 


where po is the local DM density, f{v) is the DM speed distribution normalised to unity and 
velocities are expressed in the detector reference frame. The integration over the DM velocity 
is computed from the minimum DM speed needed to induce a recoil of energy Er, Vmin = 
, where pat is the WIMP-nucleus reduced mass, to the escape velocity 
Vesc also in the detector reference frame, above which DM particles are not gravitationally 
bounded to the Galaxy. The DM-nucleus differential cross section do/dER is computed from 
the Lagrangian that describes the interaction of a given WIMP with ordinary matter and 
encodes the microscopic interactions of DM particles and quarks inside nucleons. 

® Typically given in units of counts/day/kg/keV. 
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As mentioned before, direct detection experiments usually assume the SHM, which de¬ 
scribes the Galaxy as an isotropic isothermal sphere. This model provides a useful and 
common framework to compare different experiments, but definitely it is not a realistic de¬ 
scription of our Galaxy. Due to this fact, we will use the mean speed distributions of Ref. [38] 
derived separately for three possible choices of the DM density profile, namely a Navarro- 
Frenk-White (NFW) profile, an Einasto one and a Burkert profile®. These distributions are 
extracted from self-consistent DM phase space distributions that are in agreement with the 
astrophysical observables of the Galaxy and thereby, we expect them to give a more realistic 
description of the speed distribution of DM particles in the Solar neighbourhood. For com¬ 
pleteness and to make it easy for the reader to compare with other results in the literature, 
including those provided by the experimental collaborations, we will use the SHM as well. 


In general, the DM-nucleus scattering cross section can be split into a SI and a SD 
contribution. The total cross section can be then calculated by adding these terms coherently, 
using nuclear wave functions. Thus, the differential cross section is given by 
da TTijv 


dER 


u, 


SI,N j^2 _i_ SD,N t^2 


SD 


(Er)) , 


( 2 . 2 ) 


where Uq ’ and Ug ’ are the SI and SD components of the cross sections at zero momentum 
transfer, and the form factors F^j sj;){Er) account for the coherence loss which leads to a 
suppression in the event rate for heavy nuclei in the SI and SD contributions (see Ref. [32] 
for a complete description of these prescriptions). Thanks to the Fermi’s Golden Rule in 
the Born approximation, we can factorize out all the energy (momentum) dependence of the 
scattering of DM off a nucleus inside a form factor (g). This allows us to express the zero 
momentum transfer cross sections as [32,33] 

2 


cr. 


S1,N 

'O 


< 7 , 


SD,N 

'O 


^Ep 


Z + {A-Z)^-f- 

Jp- 


a, 


SI 


+ 5.^ 

CLp _ 


1 2 


4 J+ 1 
3 J 


a. 


SD 


(2.3) 


where fip is the proton-WIMP reduced mass; Sp and Sn are the expectation values of the 
total proton and neutron spin operators; fp, fn and Up, an are the effective WIMP couplings 
to protons and neutrons in the SI and SD case, respectively^. The target material is defined 
by the following parameters, the atomic number Z, the mass number A, and the total nuclear 
spin J. In these expressions, we have defined the WIMP-proton cross sections {cip^,ap^) as. 


cr, 


SI 


_ 4 2 


2 

‘pJp 


®The mean values of the escape velocity, Sun’s velocity, and local DM density for each profile can be found 
in tables 2, 3 and 4 of Ref. [38]. 

^Notice that for simplicity we have not included here a possible vector coupling (corresponding to non- 
Majorana DM particles), which would lead to an extra contribution in the expression for 
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cr. 


SD 


24G|n 2 2 


(2.4) 


where Gi? is the Fermi coupling constant and the spin of the DM particles is considered equal 
to that of the proton. Note that Eq. (2.3) reduces to these expressions when the quantities 
referred to the nucleus are substituted by those of the proton, i.e. using = 1/2, = 0, 

J = 1/2, A = Z = 1. 


2.1 Differential rate in SuperCDMS 


The actual energy observed by the SuperCDMS experiment is not the recoil energy presented 
in Eq.(2.1), but the total phonon energy produced by a collision event. To take into account 
such a change, we follow the prescription given in Ref. [34]. The phonon energy, Ep, can be 
expressed as the sum of the recoil energy and the Neganov-Luke phonon energy, produced 
proportionally to the charge energy, 

Ep = Er + ^eAV. (2.5) 

where AV = 41/ is the bias voltage applied to the detectors, e is the electron charge, and 
hnally Eq is the mean charge for nuclear recoils. The latter is calibrated as a function of Ep 
and then parametrised as follows, 

Eq = f{Ep) = ai + a^Ep + lO^^erf , (2.6) 


where the parameters ai depend on the corresponding detector and can be found in Ref. [34]. 
Einally, to convert the differential rate terms of the nuclear recoil energy, given in Eq.(2.1), 
to the differential rate in total phonon energy, we apply the following change of variables. 


dR 

dEp 


(Ep) 


dEji 

dEp 


(Ep) 


dR dEff, ^ . 

{Er{Ep)) X {Ep), 


dER 


1 - fiEp) X 


AE 

sF 


(2.7) 


In figure 1, we show the differential rate in a Ge experiment as a function of the recoil 
energy (left panel) and the total phonon energy (right panel) for a SI cross section off nucleons 
cjp jj = 10“® pb, neglecting contributions from the SD interactions. In the right panel, we have 
used the charge model of the T2Z1 detector as it is given in Ref. [34]. We have generated the 
differential rates for distinct DM masses and different speed distributions to show explicitly 
the impact that the velocity distribution choice has in the differential rate. Our SHM results 
(solid lines) are in very good agreement with that of the SuperCDMS Collaboration when 
comparing these results with figure 2 of Ref. [34] . As it is shown, generally, the variation of 
the speed distribution tends to predict higher rates in the low mass region (all DM masses 
used in this figure would fall in such category) with respect to the SHM case. This would 
induce more stringent bounds as we will see later on. 
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Figure 1: Differential spectrum for DM-nucleon scattering on Ge in terms of the recoil energy (left) 
and the total phonon energy for the T2Z1 charge model (right) for different DM masses: = 5 GeV 

(red), = 7 GeV (blue), = 10 GeV (green) and = 15 GeV (magenta). The SI cross 
section for protons and neutrons is = 10“® pb, while the SD cross sections are set to zero. 
Different assumptions about the DM speed distribution are used: SHM (solid) and self-consistent 
speed distributions corresponding to a NFW (dashed-dotted), Einasto (dotted) and Burkert DM 
profiles (dashed) from Ref. [38]. 

2.2 Differential rate in LUX 

We address now the LUX experiment. In general, the observed energy in liquid Xe-based 
detectors is the number of photoelectrons (PE) S'!. To simulate the signal in LUX, we have 
followed the description given in Ref. [35], where the number of expected PE, v{Er), reads 

= EjiCeff{Eji)Q^, (2.8) 

where = 0.14 is the photon detection efficiency for LUX [2]. Cf,ff{Epi) is the absolute 
scintillation yield, which we have extracted from Ref. [36]. The signal rate in number of 
photoelectrons n is then given by 

TO poo TO 

^ = / dEn^ X Poiss {n\u{ER)) , (2.9) 

an Jq dEji 

where Poiss (nji^(Ej:j)) is a Poisson distribution with expectation value i'(Er), and the dif¬ 
ferential rate in Er is given by Eq. (2.1). Taking into account the finite average single¬ 
photoelectron resolution (Jpmt = 0.37 PE of the photomultipliers [37], the resulting ^l- 
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Figure 2: Differential spectrum for WIMP-nucleon scattering on Xe in terms of the recoil energy 
(left) and the scintillation signal S'! in PE (right) for different WIMP masses: = 10 GeV (red), 

= 20 GeV (blue), = 50 GeV (green) and = 100 GeV (magenta). The cross sections 
and speed distributions used are the same as in Figure 1. In the right panel, we have included the 
efficiency of the LUX detector. 

spectrum is given by 

JO °° JO 

= ^Gauss(51|n, V^no-PHT) ^ C(5'l), (2-10) 

n=l 

where C('51) is the acceptance of the applied cuts. Gauss(S'l|n,-y/ncrpMT) denotes a normal 
distribution with mean n and standard deviation ^/na■pMT■ 

In figure 2, we show a set of differential rates for different DM masses and speed distri¬ 
butions, as a function of the recoil energy (left panel) and SI (right panel) in the DM search 
energy window. The non-monotonic behaviour of the differential rates in the right panel is 
due to the efficiency of the LUX detector. As we can note, now the spectra are flatter than 
for Ge (see figure 1) as a consequence of the heavier Xe isotopes. The same tendency as 
before is observed, the predicted rates are higher for the self-consistent speed distributions 
than that predicted for the SHM. In this case, since the DM masses used here are heavier 
in some cases than that considered for SuperCDMS, we see that for a DM of 50 GeV the 
effect of the speed distribution starts to be negligible. This effect will also appear in the LUX 
upper limits. 
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3 Upper limits for SuperCDMS and LUX 


To extract the upper limits, we use the Yellin’s maximum gap method [39] at 90% C.L., which 
has been shown to produce good results for both LUX [40-48] and SuperCDMS [45-48] . This 
statistical tool allows to extract the upper bounds in presence of an unknown background in 
terms of the expected DM events in the signal region. To such end, we calculate the total 
number of expected DM events using the procedure outlined in Appendix A for different 
speed distributions as well as SD structure functions. Notice that this method does not 
permit to combine the information from both experiments at the same time, however, there 
has been other proposed methods, like the ones presented in Refs. [25,48], that allow such a 
combination. In the following, we provide the information about the SuperCDMS and LUX 
experimental setups used to calculate the tabulated data attached to this work, and therefore, 
the upper limits. 

• SuperCDMS: The SuperCDMS Collaboration has reported a first search for DM using the 
latest background rejection capabilities for an exposure of 577 kg-days [1]. This analysis relies 
on a low-threshold performance of the detector, and hence, it was carried out only for DM 
masses below 30 GeV. Eleven events were observed after unblinding, allowing to set the most 
stringent upper bound on the SI DM-nucleon cross section for masses approximately below 
6 GeV. For this analysis, only seven detectors (TlZl, T2Z1, T2Z2, T4Z2, T4Z3, T5Z2 and 
T5Z3) had significant sensitivity to very low energy recoils, and then these were employed to 
set the upper limit. The remaining detectors were used to veto events that scatter in multiple 
detectors instead, since multiple scatter event topology has a negligible probability to occur 
for DM. 

We calculate the expected signal for each detector separately using the efficiencies, expo¬ 
sures and charge models specified in Ref. [34] in terms of the total phonon energy given in 
Eq. (2.7). We also include an estimation of the energy resolution, considering the 1.3 keV 
activation line at a total phonon energy of approximately 3 keV from Fig. 3 of Ref. [1], 
to be energy independent and AEp = 0.25 keV. To find the upper limits using the maxi¬ 
mum gap method we have taken into account the eleven observed events in the signal region 
Ep =[2-12.5] keV. 

• LUX: The LUX Collaboration reported null results in the search for DM particles for an 
exposure of 10065 kg days, which allowed them to place the most stringent upper limits on the 
SI interactions of DM off nucleons above 6 GeV approximately. The signal region was set to 
be 2-30 PE in 51. The acceptance of the detector is shown in the bottom of Fig. 1 of Ref. [2] 
and we add an extra 1/2 factor to account for the 50% of nuclear recoil acceptance. We use 
the 51 single PE resolution apuT = 0.37 PE [37], a 14% of photon detection efficiency, and 
the absolute scintillation efficiency digitized from Ref. [2]. 
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Figure 3: Upper bounds on the SI proton cross section as a function of the DM mass. Black and 
grey curves are the official upper bounds from the LUX and SuperCDMS Collaborations, respectively. 
Blue and green curves are our calculations of the 90% C.L. upper limits for the LUX and SuperCDMS 
experiments respectively using the Yellin’s maximum gap method. 

In figure 3, we show the upper bounds on the SI proton cross section versus the DM 
mass for SuperCDMS and LUX experiments. In black (grey), we have plotted the official 
limit of the LUX [2] (SuperCDMS [1]) Collaboration, whereas in blue and green, respectively, 
we represent our findings of the upper limits at 90% C.L. using the Yellin’s maximum gap 
method. For consistency, our results have been obtained using the same speed distributions 
as in the case of the collaborations. In the case of LUX, to extract the upper bounds we 
have considered zero candidate events observed, even though the results of the collaboration 
showed one candidate event marginally close to the background region in the log^o(‘S'2/51) — 
51 plane. This fact is reflected as a slight disagreement between our results and those 
from the collaboration for masses below 30 GeV. On the contrary, for heavier masses, the 
agreement between our findings and the official limit is excellent. For SuperCDMS, the small 
discrepancies between our results and those of the collaboration are only due to the statistical 
method applied, since the latter used the Yellin’s optimum interval method [39]. 

3.1 Limits on the SI cross section 

Let us start with the results for the SI interactions of DM off nucleons. In figure 4, we show 
the upper limits for the LUX (blue) and SuperCDMS (green) experiments. All these curves 
have been derived under the usual assumptions, fn/fp = 1 and = 0. Different type 
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Figure 4: Upper bounds on the SI proton cross section as a function of the DM mass for the LUX 
(blue curves) and SuperCDMS (green curves) experiments. The DM speed distribution used in each 
case is: SHM (solid), NFW (dashed-dotted), Einasto (dotted) and Burkert (dashed). Darker colours 
correspond to po = 0-3 GeV/cm^ while lighter colours to the mean value of the local density extracted 
from Ref. [38]. 

of curves correspond to different choices of the DM speed distribution. For the SHM (solid 
curve), we have used the same values of the speed distribution as in Refs. [1,2] with a local 
DM density po = 0.3 GeV/cm^. On the other hand, for the remaining profiles, our results 
have been obtained with the speed distributions from Ref. [38] and are shown as dashed- 
dotted curves for the NFW profile, dotted for the Einasto profile, and dashed for the Burkert 
profile. In addition, we have rescaled all the upper bounds according to different values of 
the local DM density. We depict in darker colours the upper limits in which we have used 
the same value of the local density as in the SHM case, in order to show how the shape of 
the DM velocity distribution impacts on the upper bounds. In contrast, curves with lighter 
colours correspond to upper limits derived using the mean values that the authors of Ref. [38] 
obtained for each of the above mentioned profiles, namely, = 0.41, = 0.42 and 

= 0.41 GeV/cm^. 

As expected, the impact of the assumed DM density profile on the upper bounds is 
specially notable in the low mass region, where small changes in the tail of the speed dis¬ 
tribution function can enhance or decrease prominently the number of expected DM events 
above the threshold [11-16]. Therefore, in the < 30 GeV region, our upper limits using 
well-motivated DM profiles and a self-consistent extraction of the speed distribution tend to 
be more stringent than using the SHM, as expected from figures 1 and 2. These differences 
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Figure 5: Upper bounds on the SD proton (right) and neutron (left) cross sections as a function of 
the DM mass for the LUX (blue curves) and SuperCDMS (green curves) experiments. The DM speed 
distributions used in each case are: SHM (solid), NFW (dashed-dotted). Darker colours correspond 
to po = 0.3 GeV/cm^ while lighter colours to the mean value of the local density extracted from 
Ref. [38]. The line width of each curve denotes the use of different nuclear structure functions. For 
LUX, the thick curve is chiral EFT, the medium is the Bonn A potential, and the thin is the Nijmegen 
potential. For SuperCDMS, the thick is the chiral EFT, the medium is the Dimitrov et al. calculation 
and the thin is the Ressell et al. calculation. See text for more details. 

can be as high as one order of magnitude for masses below 10 GeV. Notice that the mean 
values of Vesc we have used are in the three cases below 544 km/s, the value used for the 
SHM. This means that the difference in the upper limits cannot be ascribed to the different 
escape velocities assumed. In fact, this shows that the slope of the speed distribution at 
high velocities is smaller in the case of the NFW, Einasto and Burkert profiles. Needless to 
say that the values used for these distributions are subject to important uncertainties which 
actually in some cases would push these upper bounds closer to the SHM result. 

3.2 Limits on the SD cross section 

Now, we move to the case of SD interactions. Our results for the upper limits on the SD cross 
section of DM off protons (left panel) and neutrons (right panel) for the LUX and SuperCDMS 
experiments are displayed in figure 5. We show the SHM (solid) and NFW (dashed-dotted 
curves) cases, where we have included different SD structure functions. Namely, we have 
employed the results obtained with chiral effective field theory (EFT) currents [49] for Ge 
and Xe isotopes (thicker lines), we have also used the Bonn A potential [50] for Xe and the 
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Figure 6: LUX combined SI and SD upper limits using a NFW profile for different masses, and 
o-nicip = 1 (left) and an/ap = —1 (right). The colour convention is as in Fig. 2. Horizontal and 
vertical lines represent the limits for SI and SD interactions only, respectively. 

Dimitrov et al. calculation [51] for Ge (medium lines), and finally, the Nijmegen potential [50] 
for Xe and the Ressell et al. calculation [52] for Ge (thinner lines). The importance of this 
kind of uncertainties on direct detection experimental results has been shown to be important 
in Ref. [9], which motivates us to include the upper bounds for the different calculations of 
the SD structure functions in the literature. We have checked that these results are in good 
concordance with those of Ref. [53] which were performed using a profile likelihood analysis. 

We note from figure 5 that due to the nuclear spin of both Ge and Xe, the constraints 
from the neutron component of the cross section is substantially more stringent than that 
of protons. Nonetheless, the difference between distinct calculations of the SD structure 
functions is smaller in the neutron case. These two facts together imply that the exclusion 
of models from LUX and SuperGDMS results is not going to undergo strong variations from 
the choice of the SD structure functions. By contrast, from the left panel of figure 5 one 
can see that for both, LUX and SuperGDMS, the change in the owing to different SD 
structure functions can be of the order 10. Of course, this might be particularly important 
for models in which the DM particle couples to protons rather than to neutrons. 

3.3 Combined limits for SI and SD interactions 

In this section, we derive the LUX limits for a specific choice of the DM mass, and then 
we calculate the upper bounds on a^p , and the combined limits. We have assumed 
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fn/fp = 1 and anjap = ±1. We perform this exercise to evidence the difference of taking 
into account both cross sections at the same time instead of considering each one separately. 

In Figure 6, we show the LUX combined bounds for a NFW profile, using pQ = 0.3 GeV/cm^, 
and for the aforementioned neutron to proton ratios for the SI and SD interactions. The struc¬ 
ture functions used in the SD case are those from the chiral EFT [49]. Besides, we have fixed 
the DM mass to: 10 GeV (red), 50 GeV (green) and 100 GeV (magenta curves). Horizontal 
lines are the SI limits considering zero contributions from the SD interactions, while vertical 
lines represent the opposite, i.e. SD limits with no contribution from the SI interactions. 
The rectangles delimited by these two kind of lines represent the cases in which the limits 
are considered for each interaction separately. Gonversely, when one takes into account both 
interactions at the same time, each rectangle bends and a new excluded-to-be region enclosed 
by the rectangles and the new limits, labeled by SI -|- SD, arises. Of course, a source of vari¬ 
ation of these combined limits is given by the change of the ratios fn/fp for SI interactions, 
and an/dp for the SD ones. For the latter, in Figure 6 we have plotted two different values, 

= 1 in the left panel, and a^/op = — 1 in the right panel. The difference among these 
two results highlights the importance that these quantities have in the upper bounds. Re¬ 
garding the SI interactions, the effect can be also remarkable, specially when fn/fp gets close 
to —0.7, i.e. the Xe-phobic region, but also for |/n//p| 1- 

Indeed, the aforesaid excluded regions are not taken into account if one applies the LUX 
results for the SI and SD interactions separately. Obviously, to extract information to rule 
out these regions it is necessary to compute the combined upper limits, but since the SI 
contribution is normally dominant, experimental collaborations only show results for one 
cross-section at a time. However, our results pose a question about the importance of the 
inclusion of such information to probe consistently DM models. As already stated, usually the 
SI contribution to direct detection experiments is dominant over that of the SD, because the 
former is summed coherently for all nucleons, hence scaling as (see Eq. 2.3). Nevertheless, 
there are scenarios in which the SD contribution can be of the same order than that of the 
SI, or even dominate over it [54]. In the following, we will show that proper combined limits 
are an important ingredient to take into account when analysing models. To such end, we 
will study two well motivated and popular DM models, the NMSSM and a Z' portal. 


4 Next-to-Minimal Supersymmetric Standard Model 

The NMSSM is one of the most popular supersymmetric models (see for instance Ref. [55] 
for a review and Ref. [56] for a DM related analysis). In this model, the DM candidate is 
the lightest neutralino (x?) which is, in general, a mixture of the superpartners of the gauge 
bosons and the Standard Model Higgs, plus a singlet scalar which is added in order to solve 
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Parameter 

Scan 1 

Scan 2 

Scan 3 

Ml 

1 - 200 

1-40 

1 - 200 

M 2 

200 - 1000 

200 - 1000 

700 - 1000 

tan/3 

4-20 

4-20 

2-50 

A 

0.1 - 0.6 

0.1 - 0.6 

0.001 - 0.1 

K 

0 - 0.1 

0 - 0.1 

0.1 - 0.6 

A\ 

500 - 5000 

500 - 5000 

500 - 1100 

An 

-50 - 50 

-30-0 

-50 - 50 


110-250 

160 - 250 

200 - 400 

Ms 

1500 

1500 

1500 

Au 

3700 

3700 

3700 

Ad 

2000 

2000 

2000 

Ae 

-1000 

-1000 

-1000 

ruf = 

300 

300 

300 

^Q, = 

1500 

1500 

1500 


Table 1: Input parameters for the series of NMSSM scans. Masses and trilinear parameters 
are given in GeV. All parameters are defined at the EW scale. 


the /x-problem. Since neutralinos are Majorana particles they contribute to both, SI and SD 
interactions, and thereby are perfect candidates to search for the aforementioned cases. 

Due to the complexity of the model, we will only focus on light neutralino masses, namely 
below m^o < 200 GeV. Furthermore, the speed distribution impact on the upper limits 
is maximised in this region. Several analyses have addressed light neutralino DM in the 
NMSSM [57-70]. 

To look for phenomenologically viable solutions of this model, we have carried out a 
series of scans over the parameter space with MultiNest 3.9 [71-73]. To this aim, we have 
defined a likelihood function whose inputs are the neutralino relic density, the SM Higgs mass, 
BR(Hs —)• and BR(6 —>■ sj), which are taken as gaussian probability distribution 

functions around the measured values with 2a deviations. The gluino soft mass is fixed 
to M 3 = 1500 GeV, the trilinear parameters to Au = 3700 GeV, Ad = 2000 GeV, and 
Ae = —1000 GeV, and the soft scalar masses of sleptons and squarks to = 300 

GeV and rUg, = nifj, = = 1500 GeV, respectively, where the index i runs over the 

three families. This conservative choice of the squark masses is motivated by the LHC null 
searches. Also note that despite the high trilinear term Ajj, the instability against charge- 
and/or color-breaking minima is avoided since the squark soft masses are at the TeV scale [74]. 

We have performed three scans over the parameter space where the different input pa- 
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rameters are varied according to Tablet. The neutralino relic abundance has been calcu¬ 
lated with micrOMEGAs 3.6.9 [75]. We have used NMSSMTools 4.1.2 [76-78] to compute 
the NMSSM mass spectrum, the masses of the Higgs bosons including full two-loop con¬ 
tributions, and the relevant low-energy phenomenology observables. LHC measurements 
of the Higgs properties are included by constructing the distribution from Ref. [79] 

and allowing 3a deviations. Limits on the velocity averaged cross section of DM parti¬ 

cles from the gamma-ray null observation of dwarf spheroidal galaxies (dSphs) are also 
taken into account, following the procedure sketched in Ref. [80] for the Pass 8 data [81]. 
In addition, we have implemented the recent measurements of the branching ratio of the 
Bg —>• process by the LHCb [82] and CMS [83] Collaborations, which collectively yield 

1.5 X 10“® < BR(Rs —)■ < 4.3 x 10“® at 95% CL. For the 6 —)• sy decay, we have 

considered the 2a range 2.89 x 10“'^ < BR(6 —)• sy) < 4.21 x 10“^, which takes into ac¬ 

count theoretical and experimental uncertainties added in quadrature [84-88]. We have also 
imposed 0.85 x 10“'^ < BR(R''" —> < 2.89 x 10“^ [89]. 

Finally, we have set an upper bound on the neutralino relic abundance, < 0.13, 

consistent with the latest Planck results [90]. Besides, we have considered the possibility that 
neutralinos might only contribute to a fraction of the total relic density. To deal with these 
cases, the fractional density, ^ = min[l, D^o/i^/0.11], have been introduced to account for the 
reduction in the rates for direct and indirect searches®. 

4.1 Combined limits for different speed distributions 

In this section, we will discuss how the inclusion of direct detection upper limits taking into 
account the SI and SD interactions at the same time can have an impact on the viability of 
DM scenarios in the NMSSM that otherwise will be allowed. More precisely, we will use the 
results from previous sections to illustrate how different ways of including direct detection 
bounds can affect the parameter space of the NMSSM. We will also show the impact of 
the astrophysical uncertainties from the speed distributions, leaving to the next section the 
analysis of the effect of these uncertainties on the SD structure functions. 

To show schematically how different solutions are excluded depending on the assumptions 
considered when evaluating the limits, we will use the following procedure: 

(1) SI: We will apply first the SuperCDMS and LUX bounds only for the SI interactions and 
f/fp = l to the viable solutions found, using the SHM. This bound correspond to the usual 
bound published by the SuperCDMS and LUX Collaborations. 

(2) SI/SD: Then, we will add to the previous limits the bounds that only take into account 
the SD interactions, also for the SHM. This step is equivalent to apply the combined limits 

^Although the introduction of ^ for certain cases of multicomponent DM, such as self-interacting DM, 
might not be consistent. 
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Figure 7: Excluded points for neutralino DM in the NMSSM in the plane using different 

assumptions on the LUX and SuperCDMS limits and different speed distributions. The blue (green) 
curve represents the LUX (SuperCDMS) bound only for SI interactions and for the SHM. 
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generating the rectangles shown in Fig. 6. 

(3) SI + SD: Next, instead of the previous considerations we will calculate the bounds taking 
into account all the contributions, SI and SD and the specific neutron to proton ratios. This 
is also done for the SHM. 

(4) Profile: Same as (3) but for the NFW, Einasto and Burkert velocity distributions. 


Notice that with the purpose of isolating the effect that speed distributions have on the 
exclusion regions, in none of the previous steps we have varied the value of the local density, 
Pq. In fact, if we change the local density value in agreement with those extracted in Ref. [38], 
according to Fig. 4, the upper bounds become more stringent. 

The set of plots shown in Figure 7 depicts how the previous steps affect the exclusion 
regions of neutralinos in the NMSSM. All the plots show the points excluded by the Super- 
CDMS and LUX results in the plane . In the first row left-hand side plot, we 

display the solutions excluded when the SI proton cross section weighted by ^ is above the 
nominal SuperCDMS and LUX limits. This corresponds to step (1), as labelled in the plot. 
In the same row, the right-hand side plot, includes also the upper bounds both for SD proton 
and neutron interactions. As it can be seen, when compared with respect to the left-hand side 
plot, a new population of excluded points has appeared around m^o ~ 30 GeV. These points 

A1 

correspond to solutions in which the SD neutron cross section is above the LUX bound®. 

Next, we move to the middle row of Figure 7. In the left panel, we note that many 
solutions previously allowed by LUX are now ruled out. These points fall precisely in the 
regions labelled by SI -|- SD in Fig. 6. All the points that have appeared at this step have 
a similar contribution from the SI and SD interactions to the expected number of events in 
LUX. Therefore, all these solutions require a careful calculation of the LUX upper bound 
to be correctly excluded. This nicely highlight the importance of using combined limits and 
showcases how the upper limits should be compared point per point in the parameter space. 
This task can be achieved very quickly with the help of the tables provided in this work. 

The remaining plots correspond to step (4). As previously discussed, the impact of con¬ 
sidering different velocity distributions is notable for neutralino masses below approximately 
50 GeV, as expected from the results shown in Fig. 4. This light mass region is actually 
one of the favoured regions that explain the Fermi-GeV excess in terms of annihilating DM 
particles [91-100]. Moreover, this excess prefers a NFW-like profile; hence, the use of realistic 
speed distributions according to the Galactic DM profile is required to suitably constraint 
possible DM explanations of this excess with direct detection upper bounds [103]. 

Thus, our findings point out that in order to explore robustly the parameter space of the 

^Remind that Xe-based experiments are much more sensitive to the neutron contribution of the SD inter¬ 
actions since the total spin of these nuclei are dominated by neutrons. 
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m^o (GeV) m^o (GeV) (GeV) 


Figure 8: Excluded points for neutralino DM in the NMSSM in the plane ^(T?o_p-m^o using different 
speed distributions and SD structure functions. The blue (green) curve represents the LUX (Super- 
CDMS) bound only for SI interactions and for the SHM. Light blue points denote solutions in which 
the contribution to the expected number of events in LUX is dominated by SD interactions. 
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NMSSM, one has to include properly the SuperCDMS and LUX limits. It is worth noting 
that steps (1) and (2) assume a fixed neutron to proton ratio for both interactions, which 
does not coincide exactly with those of neutralinos, when we apply step (3) this subtlety is 
considered using the exact values that neutralinos have for these quantities. Therefore, these 
ratios can affect the exclusion of some of the previously allowed points. 

4.2 Combined limits for different SD structure functions 

As seen before, many solutions found in our scans predict a SD contribution similar to that 
of the SI. In Section 3.2, we have shown that the use of distinct structure functions lead to 
differences in the LUX upper bounds that can be significant. Therefore, we must expect that 
the choice of the structure functions is going to impact on the excluded regions of neutralinos 
in the NMSSM. 

In Figure 8, we depict the excluded points assuming different SD structure functions 
(columns) and speed distributions (rows). For the sake of completeness, we have plotted in 
light blue those points in which the contribution to the DM expected events in LUX from the 
SD interactions is dominant. As we can observe, in Figure 8 from left to right panels, light 
blue points are generally affected by the choice of the structure functions and, in agreement 
with Figure 5, the chiral EFT structure function is able to rule out smaller regions of the 
parameter space. It is noteworthy that the Bonn A potential model of the SD structure 
functions provides the most stringent limits of Xe-based experiments, highlighted as a higher 
population of light blue points in the middle column of Figure 8 with respect to the other 
two (compare also steps (2) and (3) in Figure 5). Notice that the SD dominant solutions we 
have found accumulate around neutralino masses of 30 GeV, while those with similar SI and 
SD contributions in LUX range from 20 to 100 GeV approximately (see also in Fig. 7 the 
change from step (2) to step (3)). This means that the choice of the SD structure function 
has a notable effect on the allowed (excluded) solutions for NMSSM neutralinos below 100 
GeV, which is precisely the region where most of the explanations of the Fermi-GeV excess 
in this model lie [66,68-70,104]. 


5 Z' portal 

We now turn our attention to another popular and well-motivated extension of the SM, 
the Z' portal. In this scenario, a hidden sector communicates with the SM via a gauge 
boson, provided that the SM is enlarged with an extra abelian gauge group [105]. The 
phenomenology of these constructions is very rich, and ranges from colliders to direct and 
indirect searches for DM [106-118]. 
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Parameter 

Range 

a 

10 "^ - 1 

l^l 

10 "^ - 1 

c 

10 -^ - 1 

\d\ 

10 "^ - 1 

\h\ 

10 -^ - 1 


1 - 200 

mz' 

500 - 8000 


Table 2: Input parameters for the scan of the Z' portal model. Masses are given in GeV. 


We will consider the same construction as in Ref. [116], which is based on a certain type 
II string compactifications with intersecting branes. The key point of these scenarios is that 
gauge bosons acquire a mass through the Stiickelberg mechanism [119,120], except of course 
the one corresponding to the hypercharge, which remains massless in the phase of unbroken 
electroweak symmetry. 

We will assume throughout this section that the DM in this scenario is a Dirac fermion, 
ij}, which lives in the hidden sector. The phenomenology of this kind of DM will be driven by 
the interaction with the SM fermions through the exchange of a Z' boson. The couplings of 
the Z' to the SM fermions will be determined by specific combinations of four parameters, 
a, b, c and d [116]. We have performed a scan over the parameter space, with the parameter 
ranges given in Table 2, where h is the coupling of DM to the Z' boson. To this end, we have 
implemented the model in micrOMEGAs 3.6.9 [75]. We have taken into account that there 
exist certain combinations of the these parameters that are not allowed by the orthogonality 
and normalization conditions imposed on the Z' mass eigenvectors [116]. 

In general, ?7(1) extensions of the SM can be constrained using resonances at colliders, 
since the coupling of the Z' boson to leptons and quarks contributes to the appearance of 
dimuon, dielectron and dijet resonances. Here, we have applied the ATLAS results for high 
mass resonances decaying into a ii'^ or an e"*“e“ pair at a center of mass energy y/s = 
8 TeV and luminosities of 20.5 fb“^ and 20.3 fb“^ for dimuons and dielectrons resonances, 
respectively [121]. Furthermore, we have used the searches for dijet resonances and monojets 
plus missing energy that receive additional contributions from the presence of a Z' boson, 
both at the LHC and Tevatron colliders [122-124]. In order to implement these constraints, 
we have followed the procedure of Ref. [116]. 

It is worth mentioning that we do not impose any bound on the relic abundance of ip. 
This is a consequence of the potential complexity of the gauge structure and matter content 
of the hidden sector, which might induce other contributions that account for the correct 
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Figure 9: Excluded points for Dirac DM in the Z' portal model in the plane using different 

assumptions on the LUX and SuperCDMS bounds and different speed distributions. The blue (green) 
curve represents the LUX (SuperCDMS) bound only for SI interactions and for the SHM. 


22 









abundance such annihilations of ijj into hidden sector matter. Of course, this will affect the 
predictions of the model for indirect dark matter searches. To deal with this, we have defined 
^ as in the NMSSM, but in this case if > 0.13 then we set ^ = 1.^® 

Finally, we have imposed the upper bounds on the annihilation cross section of if: from 
dSph galaxies using the Pass 8 data of the Fermi-LAT satellite [81]. These limits are presented 
for specific SM final states, while in this model each point of the parameter space entails DM 
annihilation into a combination of different SM particles. This circumstance prevents us 
to use these Fermi-LAT results directly, but we have to weight each upper limit by the 
corresponding annihilation percentage instead. 

5.1 Combined limits for different speed distributions 

As in the NMSSM case, presented in Section 4, in this section we will evaluate the impact 
of the speed distribution on the experimentally allowed solutions of the Z' portal with Dirac 
DM. We will follow the same steps described before for the NMSSM. 

In figure 9, we show the results for the Z' portal model using the same steps as in figure 7. 
Notably, from step (2) to (3) many new excluded solutions appear below the LUX SHM limit 
(blue curve). Unlike in the NMSSM, all these ruled out points correspond to solutions in which 
the ratio fn/fp is different from 1, and in many of these cases, the neutron contribution to 
the SI cross section is much higher than that of the proton. For this reason, many points that 
lie far below the LUX bound on the Sl-proton cross section are excluded. In addition, some 
of the solutions that were excluded in step (2) are now allowed. These solutions correspond 
to those with fn/fp ~ —0.7, the region in which the sensitivity of Xe-based experiments 
substantially decreases, the so-called Xe-phobic region. Also, solutions with |/n//p| < 1 are 
now allowed, since the neutron contribution does not reach the size of the proton one, and 
hence, the upper bound weakens. 

Regarding the impact of the velocity distribution, one can see in figure 9 that from step 
(3) on, more points are excluded in the low mass region. Although in this case, the difference 
is less visible than in the NMSSM, because the density of points for masses below 30 GeV is 
not very high. Furthermore, the fact that in this model fn/fp can virtually take any value, 
makes more difficult to identify points where the effect of the speed distribution is remarkable. 
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Bonn A potential, Dimitrov etal. 


Nijmegen potential, Resselt et at. 





(GeV) 


Figure 10: Excluded points for Dirac DM in the Z' portal model in the plane using 

different speed distributions and SD structure functions. The blue (green) curve represents the LUX 
(SuperCDMS) bound only for SI interactions and for the SHM. Light blue points denote solutions in 
which the contribution to the expected number of events in LUX is dominated by SD interactions. 
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5.2 Combined limits for different SD structure functions 


In figure 10, we show the results for the Z' portal model using the same speed distributions 
and SD structure functions as in figure 8. Notice that the absence of light blue points (like 
those appearing in figure 8 for neutralino DM) means that in this case we have not found any 
point of the parameter space in which the SD contribution to the total number of expected 
events in LUX dominates. Consequently, we observe that the change in the SD structure 
functions employed to impose the limits does not affect the exclusion of the solutions found. 


6 Conclusions 

In this paper, we have performed a thorough study of the impact of the assumptions usually 
made to compute DM direct detection upper bounds. We have used the NMSSM and a 
Z' portal DM models to illustrate how the use of direct detection limits, including SI and 
SD interactions at the same time, different neutron to proton ratios, distinct DM density 
profiles as well as different SD structure functions, has important implications on the excluded 
parameter spaces. 

We have extracted the LUX and SuperCDMS upper bounds under many different as¬ 
sumptions about the speed distribution and SD structure functions. For the former, we have 
used the velocity distributions derived in Ref. [38] for three different choices of the DM den¬ 
sity profile. In Ref. [38], the authors computed the speed distribution in a self-consistent way 
from the gravitational potential of the Milky Way. We have employed these distributions, to 
compute upper limits that rely on a realistic description of the DM halo, and hence, produce 
results consistent with a wide range of astronomical observations. As expected, the differ¬ 
ences in the upper bounds from the speed distribution focus mainly on the light mass range, 
iTT'DM ^ 50 GeV. Another source of uncertainty for direct detection experiments are the SD 
structure functions. By using the different results in the literature, we have explicitly shown 
the variation in the upper limits due to the choice of the SD structure functions. 

In order to incorporate all these ingredients when analysing models, in each of the men¬ 
tioned cases, we have carried out a decomposition of the expected DM signal in both exper¬ 
iments, LUX and SuperCDMS, in terms of several expressions that can be used to derive 
upper limits using statistical tools based on the total number of expected events. In this 
paper, we have used the maximum gap method, and obtained tabulated data, that is suitable 
for this method, for each of the different combinations of these ingredients. This data set 
allows to constrain any generic model with both SI and SD interactions, for protons and 

^°This is equivalent to consider that whatever hidden processes there could be behind the relic density, 'i/i 
particles can account for the 100% of the DM in the Universe. 
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neutrons. The data and python scripts that exemplify their use can be downloaded from this 
site http://goo.gl/lCDFYi. 

As case studies to apply our findings, we have chosen the NMSSM and the Z' portal 
with Dirac DM. We have performed a series of scans over the parameter space of both 
models imposing the latest experimental limits. More specihcally, we have applied to the 
resulting points the constraints from the LHC and Tevatron, Planck, low energy observables 
and Fermi-LAT bounds from dSph galaxies. Then, we have imposed how the inclusion of 
the upper limits from LUX and SuperCDMS taking into account the SI and SD interactions, 
and different assumptions about the halo profile and SD structure functions, to analyse the 
effect of each of these factors on the allowed/excluded regions of the parameter space. Our 
findings point out that a more careful implementation of direct detection limits is needed 
when analysing models in light of the current experimental constraints. More specifically, in 
the NMSSM, the nature of the solutions found requires a combined implementation of the SI 
and SD bounds at the same time, since many solutions have similar contributions from both 
types of interaction. For this reason, in this model, not only the halo profile choice, but also 
the SD structure functions used to derive the SD contribution play an important role. In 
the case of the Z' portal, the necessity of applying combined limits arises from the neutron 
and proton contributions to the SI interactions. This model, with a Dirac DM candidate, 
shows strong deviations from the /„ = fp equality. This situation is translated into huge 
variations of direct detection upper limits that should be considered for a careful analysis of 
the parameter space. 

The amount of experimental data nowadays is enormous. Huge experimental and theo¬ 
retical efforts are being made to interpret these data in terms of many different models. With 
this purpose, theoreticians are performing very fine calculations of observables that can be 
constrained in light of the current data. In this sense, our results provide a missing piece to 
this effort, making easier the implementation of direct DM detection limits using more real¬ 
istic assumptions; namely the SI and SD contributions at the same time, speed distributions 
in agreement with astrophysical observations, and different SD structure functions. 
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A Extraction of upper limits using the maximum gap method 


In this appendix, we give detailed information about the extraction of the direct detection 
bounds for SuperCDMS and LUX using the data provided in this work. Notice that this data 
set can be used with two different, nevertheless equivalent, purposes. Its first application is to 
calculate the upper limits under a number of assumptions, in the same way the collaborations 
present their results for the SI proton cross section as a function of the mass (see for instance 
Fig. 4). The second application is to test if a given point in the parameter space of a model 
is allowed or excluded, as we have done to obtain Figs 7-10. 


As stated in the Introduction, by means of these tabulated data one can determine the 
expected number of DM events in SuperCDMS and LUX in an easy and quick way. Then, 
with these number of events, we can estimate the upper bounds on DM interactions with 
protons and/or neutrons, using the Yellin’s maximum gap method [39], which produces ex¬ 
cellent results (see Fig. 3 for a comparison between our results and those of the experimental 
collaborations). With the attached data, one can calculate the expected number of events in 
specific energy windows and the total number expected for a given experiment. The former 
and the latter are necessary to apply this method to SuperCDMS (for LUX, we assume zero 
observed candidate events and then we will only need the total number of expected events). 


Let us start considering a generic direct detection experiment, denoted by Z. It is straight¬ 
forward to show that the number of expected DM events for the SI interactions can be written 
as, 


Nsi = cr: 


SI 




(A.ll) 


where the functions enclose all the information about the DM halo, form factors, 

energy resolution, as well as the transformation from the actual measured energy to recoil 
energy (see Section 2). These SI functions have been normalized to a reference cross section 
of 10“® pb. Note that these functions have been calculated for a specific energy window, and 
then, are not valid for arbitrary energy intervals. Furthermore, we have reabsorbed all the 
constants, different for each function, namely {A — Z)^ and Z{A — Z), and the exposure 
of the detector in the definition of the functions. 


For the SD interactions, one can also perform the same exercise yielding. 


Nsd 


a, 


SD 


^ " (^) ) ’ 




Fz\rn^) 


(A.12) 
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where the superindices of the functions correspond to the isoscalar and isovector 

decomposition of the SD structure functions. These SD functions have been normalized to 
a reference cross section of 10“^ pb. Therefore, the total number of events in a generic DM 
experiment can be calculated just by summing the SI and SD contributions. Up to some 
constants, these functions are the same as those presented in Ref. [29] for Oi and O 4 in the 
framework of the effective field theory for DM interactions. 

The maximum gap method finds the most constraining gap between observed events for 
setting an upper limit. This is, the energy range between two observed events with the 
highest number of events predicted by a model. Then, to calculate the bounds we determine 
the probability of this gap being smaller than a particular number of events, x [39], 


Co{x,n) 


{kx 

k =0 


/i)^e 



^ — kx J ' 


(A.13) 


where /i is the total number of events in the whole search energy window, this is Nsi + Nsd 
from Eqs. A.11 and A.12, and m is the greatest integer < fi/x. First, note that if there is any 
observed candidate event for the experiment considered in the search window then x < fi, 
and the functions introduced in Eqs. A.11 and A.12 to calculate x and /r are not the same, 
because the energy ranges used in each case are different. Second, for zero observed events 
= X hy dehnition of the maximum gap [39], and the previous equation reduces to. 


C'o(x) = 1 - e 


(A.14) 


In the following, we will explain how to use the above mentioned expressions to find the 
upper limits and to test points in the parameter space of a given model, for the specific cases 
of SuperCDMS and LUX. 

• SuperCDMS bounds: Since SuperCDMS observed eleven candidate events in the DM 
search region [1], we have to define the functions for different energy ranges. All the functions 
have been extracted accordingly to the experimental setup of Ref. [1] (see Section 3 for more 
details) and following the prescription given in Section 2. To evaluate the limits up to 
DM masses of 30 GeV, two different mass ranges must be considered < 5.4 GeV and 
> 5.4 GeV. For each range, there is a different maximum gap and consequently we have 
to compute a distinct x, denoted as xi and X 2 - Now, to calculate fi and xi ^2 or X 2 , 
depending on the DM mass), one must use Eqs. A.11 and A.12 with the proper tabulated 
functions given in the provided datafiles (see appendix B). 

To test if a specific particle physics input characterized by /n//p, o'p^, a^jap and 
is allowed or excluded at e.g. 90% G.L. by SuperCDMS, we compute /x and xi ^2 and then 
evaluate Cq from Eq. A.13. If Cq > 0.9, the input is allowed, otherwise is excluded. On the 
other side, if we want to determine the upper limit at 90% C.L. for the SI interactions and 
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fn/fp = 1, one has to calculate ji and xi ^2 using only the SI functions and choosing fn/fp = 1- 
Afterwards, for a specific DM mass, one finds the such that Cq = 0.90. 

• LUX bounds: The calculation of LUX bounds is easier because we have considered zero 
candidate events. In this case, the calculation of Cq only requires the computation of the 
total number of expected events, /z, for the whole energy range and valid for all DM masses. 
Hence, the suitable expression for Cq is Eq. A. 14. To derive the necessary functions, we have 
used the experimental setup described in Ref. [2] (see Section 3 for more details) and followed 
the prescriptions given in Section 2. To perform the same tasks as for SuperCDMS, we follow 
the same procedure, the only difference being the dehnition of Cq. 


B Attached files 

All the datafiles obtained in this work can be downloaded from the site http: //goo. gl/lCDFYi. 
We provide separated files for the tabulated functions, introduced in Eqs. A.11 and A.12, 
corresponding to the SI and SD interactions, different speed distributions, and SD structure 
functions, for DM masses below 200 GeV. The files are named according to the following 
convention; 


SI files: <Experinient>_SI_<Prof ile> .txt, 

SD files; <Experiment>_SD_<Profile>_<FormFactor>.txt, 

where <Experiment> = SuperCDMS, LUX; <Profile> = SHM, NFW, Einasto and Burkert for 
the SHM and the speed distributions of Ref. [38]; <FormFactor> = Chiral, Ressell and 
Dimitrov for the SD structure functions of SuperCDMS from Refs. [49,51,52] respectively, 
and <FormFactor> = Chiral, BonnA and Nijm for LUX using results from Refs. [49,50]. 

For LUX, we also give results for DM masses above 200 GeV, for different SD structure 
functions, and using only the SHM, since as we have seen, different speed distributions do 
not impact on the results for masses above 50 GeV approximately. These files contain an 
extra _extended in the filename, but they must be used in the same way as the other files. 

In the header of each hie, the distribution of the columns is given. E.g. in the hie 
SCDMS_SI_SHM. txt, we hnd the tabulated functions needed to calculate /z, xi and X 2 (denoted 
by the suffixes, _Total, _MG1 and _MG2, respectively) using Eq. A.11, taking into account 
only SI interactions and for the SHM. The header of this hie reads, 

####################################################################################### 
## Functions for SuperCDMS SI interactions using a SHM halo and the Helm form factor 
## Energy intervals:MG1=[2.411, 2.764] keV,MG2=[3.561, 7.176]keV,Total=[2.0, 12.5]keV 
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## Columns:Mass,Fp_MGl,Fp_MG2,Fpn_MGl,Fpn_MG2,Fn_MGl,Fn_MG2,Fp_Total,Fpn_Total,Fn_Total 
####################################################################################### 

The LUX files only have four columns, for instance the header of the file LUX_SI_SHM.txt, 
containing all the necessary information to calculate the SI induced events for a SHM speed 
distribution reads, 

############################################################################## 

## Functions for LUX SI interactions using a SHM halo and the Helm form factor 
## Energy intervals: Total = MG = [2.000, 30.000]keV 
## Columns: Mass, Fp_Total, Fpn_Total, Fn_Total 
############################################################################## 

Finally, we also provide three example codes for LUX and SuperCDMS, namely LUXexample. py, 
SCDMSexample.py, and LUXexample2.py. These python scripts implement the expressions of 
appendix A and illustrate how to use them together with the attached datafiles. The first two 
codes evaluate if a point of the parameter space, whose specific characteristics are included in 
the scripts as an example, is allowed or excluded at 90% C.L. for each experiment separately. 

The third code reads 1000 inputs from the NMSSM parameter space (not all constraints have 
been imposed to this data), saved in testdata.dat, checks if they are excluded or not by 
LUX, and then plots the allowed and excluded points in the apL^-m^ plane. This is useful 
when one wants to compare a scan over the parameter of a given model with the combined 
LUX bounds. 
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